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include divergent trajectories. It is possible, however, to solve the stochastic 
problem exactly, but the averaging must be performed with great care. 
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I. INTRODUCTION 



The treatment of even quite simple quan- 
tum optical systems can present a significant 
technical challenge. The description of any 
isolated system can be given using a den- 
sity operator, with time evolution governed 
by the Liouville equation |[|. When the sys- 
tem of interest is not isolated, but can ex- 
change both energy and fluctuations with a 
surrounding environment, the evolution of 
the system density operator is governed by 
a master equation 0. In particular master 
equations provide a practical method to treat 
such systems but direct solution of these is 
not usually possible. It is often possible, par- 
ticularly for problems involving optical field 
modes, to map the operator master equa- 
tion onto a partial differential equation for a 
quasi-probability distribution. It may be pos- 
sible to solve this equation or to map it onto 
an equivalent stochastic process that can be 
simulated numerically. 

Mapping the quantum problem onto a 
stochastic system relies on a formal simi- 
larity between the partial differential equa- 



tion, obtained from the master equation, and 
the Fokker-Planck equation associated with 
Brownian motion. The Fokker-Planck equa- 
tion for the dynamics of a single field mode 
or harmonic oscillator is typically of the form 



d_ 
dt 



W=-— A a W-—A a *W + 
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2 da* 2 " " ' 1 dada 
where W is the quasi-probability distribution 
for the phase-space associated with the mode 
and parameterized by the complex variables 
a and a* The requirement that W be 

a real- valued function imposes the conditions 
that A* a , = A a , D a * a * = D* aa and D aa * is 
real. This equation can be mapped onto a 
pair of stochastic differential equations for 
the phase-space coordinates (also written as 
a and a*) in the form 



a = v a + 
a* =v Q * +C(t) 



(2) 
(3) 



where v a and v a * are functions of the drift 
and diffusion coefficients and Dy) ap- 
pearing in Eq. flU) and the dot denotes 
a derivative with respect to time. The 
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terms £(t) and £*(t) are stochastic fluctuating 
terms with correlation functions related to 
the diffusion coefficients. There is no unique 
stochastic representation of a given Fokker- 
Planck equation. In this paper we work with 
the Stratonovich form of the stochastic inte- 
gral M. A brief discussion of this is given in 
Appendix 

Unfortunately, not all problems of interest 
can be converted into the Fokker-Planck form 
(H) . Systems of interest in quantum nonlinear 
optics often produce equations for the evolu- 
tion of quasi-probabilities that have deriva- 
tives of higher than second order and it is not 
known how to treat these within the stochas- 
tic formalism. The usual approach is to sim- 
ply drop these terms to produce "stochas- 
tic electrodynamics". It has been shown, 
however, that this frequently used approxi- 
mation does not correctly reproduce higher- 
order correlations such as those predicted to 
occur in parametric oscillators |5j||. A sec- 
ond, more subtle, problem is that even when 
we do obtain an equation of the form ([!]), it 
might still not be possible to map this onto 
SDEs of the form (|) and (§. The difficulty 
arises when we have negative diffusion, that is 
when D aa * < \D aa \. With negative diffusion, 
we are led to SDEs in which £* cannot be the 
complex conjugate of £ and hence a* will not 
be the complex conjugate of a. It was to re- 
solve problems of this kind that the positive 
P representation was introduced 0,0-0. 

In this paper we consider a proposal by 
Yuen and Tombesi to convert the evolution 
equation for the Q quasi-probability into a 

Their idea is that the 



pair of SDEs ggfn| 
correct averages should be obtained by for- 
mal application of the Langevin method by 
simply ignoring the presence of negative dif- 
fusion. These authors applied their method 
to a single-mode evolving under the influ- 
ence of a quadratic Hamiltonian in the pres- 
ence of damping and showed that this gave 
the known evolution for this problem. In 



this paper we apply the Yuen- Tombesi ap- 
proach to the more demanding, but still an- 
alytically soluble problem of the undamped 
anharmonic oscillator fl2-14|. This model 



is known to cause difficulties with stochastic 
simulations derived from the positive P rep- 
resentation P,P|JT5[|. We find that the Yuen- 
Tombesi method gives the correct results but 
that it cannot reliably be applied to numeri- 
cal simulation of the problem. We trace the 
origin of this difficulty to the order in which 
stochastic averages and averages over the ini- 
tial phase-space distribution have to be per- 
formed. 



II. METHOD OF YUEN AND 
TOMBESI 

The method of Yuen and Tombesi was de- 
signed to deal with problems in which the 
evolution equation for the Q function dis- 
plays negative diffusion. The Q function for a 
single field mode or oscillator can be written 
in a number of forms, the simplest of which 

H 



is 



Q(a, a*) = — (a\p\a), 



(4) 



where p is the density operator for the mode. 
This distribution can be used to obtain anti- 
normal ordered moments of the annihilation 
and creation operators by integration over 
the complex a plane: 



(a n a tm ) 



j d 2 a Qa n a* m . 



(5) 



We consider systems (such as the anhar- 
monic oscillator) in which the evolution equa- 
tion for the Q function is of the form given in 
Eq. (0), with negative diffusion. This leads 
to associated SDEs in which the stochastic 
variable a*(t) is not the complex conjugate 
of a(t). As an example, consider an equation 
in which D aa * = 0. This necessarily implies 
negative diffusion associated with D aa and 
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D a * a *. We can follow the method outlined 
in Appendix ^ to obtain a pair of SDEs that 
are equivalent to our evolution equation for 

Q ET 



. , _ ld_ 



+ 



D f 



(6) 

Da* a* fa*- (7) 



It might appear that these equations are mu- 
tual complex conjugates but this is not the 
case as the two Gaussian noise terms are in- 
dependent and hence do not take complex 
conjugate values. It follows that we can- 
not interpret a and a* as mutual complex 
conjugates. The situation is reminiscent of 
that found with the positive P representa- 
tion and we employ the same notation by 
writing our stochastic variables as a(t) and 
a + (t) f7j. Anti- normal ordered expectation 
values should then correspond to stochastic 
averages of corresponding functions of a and 
a + , with a(t) replaced by a(t) and a'(t) re- 
placed by a + (t). 

We can introduce the variables a and a + 
more formally by means of the complex func- 
tion 

1 



Q(a, a 
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which is a function of a and a + but not of 
their complex conjugates. This reduces to the 
familiar Q function (|4]) when a + = a*. We 
can convert our master equation for p into 
an evolution equation for Q by making the 
substitutions 




(9) 

The resulting equation for Q will be of the 
same form as that for our Q function with a* 
replaced by a + . 



III. ANHARMONIC OSCILLATOR 



The anharmonic oscillator is one of 
the simplest analytically solvable models in 
quantum optics. The Hamiltonian for this 
model can be written in the form 



H 



uj (ala + + p(a ] a) 2 , (10) 



where uj is the natural angular frequency 
for the mode and we work with units in 
which h — 1. The term proportional to 
/i is sometimes written in normal order as 
pa^a^aa. This is the same model but with 
the uj changed to uj+p. It is convenient to re- 
move the free evolution of the mode and this 
can be achieved by working in an interaction 
picture rotating at angular frequency uj. The 
interaction picture form of the Hamiltonian 
(0) is 



Hi = p(a)a) 



(11) 



This Hamiltonian has been used in quan- 
tum optics as a model for the Kerr nonlin- 
ear refractive index. Despite its simplicity, 
it produces a number of nonclassical effects 
including squeezing JTJ]] and Schrodinger cat 
that is superpositions of coherent 
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states 

states. The accurate reproduction of these 
features, especially the cat states, presents 
a stiff challenge to a stochastic simulation 
method such as that proposed by Yuen and 
Tombesi [TJjjll|]. The fact that the model is 
analytically soluble means that we can com- 
pare the results of such simulations with ex- 
act analytical expressions. We will give an 
example of this comparison in the following 
section. In this section we present a brief re- 
view of some of the known features of the 
model. 
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It is clear from the form of the Hamilto- 
nian that it commutes with the number op- 
erator o)a. It follows that the number of ex- 
citations (or photons) in the mode will be 
conserved and that the photon number states 
\n) will be the eigenstates of our interaction 
Hamiltonian 

H^n) =n 2 n\n). (12) 

The corresponding time-evolution operator is 



U(t) = exp(-i#jt) = \n)(n\e~ 



(13) 



n=0 



and it follows that the evolution of our os- 
cillator will be periodic with period 27r//x. If 
we can expand our initial state in terms of 
the number states, then we can use this re- 
sult to solve for the time-evolved state in the 
Schrodinger picture. For example, an initial 
coherent state |a ) wm evolve to the state 



U(t)\a ) = e" |ao1 -^=e- m ^\n). (14) 
n=o Vn\ 

This state has a rich structure that can 
be seen in pictures of the associated quasi- 
probability distributions [12,14]. The state 



has a simple form at the quarter periods when 
it can be written as 



£>[7r/(2/x)]K> 
U[ic/fjL]\ao) 
U[3n/(2 l 2)}\a ) 



2 

I - «o) 
1 + i 



«o) + 



l + i, 



a ) 



in the Heisenberg interaction picture. The 
time-evolved annihilation and creation oper- 
ators are 



a(t) = U\t)aU(t) -- 
a\t) = U\t)a ] U(t) 



-ifit —i2{tta^ a 



(I 



(16) 
(17) 



where we have written the initial operators as 
a and a) . It is straightforward to use these ex- 
pressions to calculate expectation values for 
functions of a(t) and a'(t). For example, the 
expectation value of the annihilation opera- 
tor for the coherent state \cxq) is 



(a(t)Y- 



e 

-i/it 



i^t /„ I „—i2uta' &z \ _ \ 
r \OLQ\e CL\Oto) 

a exp (|a | 2 (e" j2At * - 



(18) 
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In this expression we have omitted the free- 
evolution in the form of a factor e~ luJt . This 
corresponds to working in a frame rotating 
at frequency u>, associated with our choice of 
interaction picture. All expressions in this 
paper will be given in this frame. The expec- 
tation value of a) {t) is the complex conjugate 
of Eq. (|18|) and higher order moments can 
also be calculated without difficulty. 

The evolution equation for the Q function 
can be written in the form 



12 



d n ■ 







da 



3)«*g-^-(2|a| 2 
oa 



+ 



d 2 



a* 2 Q - -^o?Q 



2 



OLq) + 



da* 2 da 2 
c*o). (15) 

Comparison with Eq. 



3)aQ 



(19) 



(P reveals that this 

The state at one quarter and three quarters equation has negative diffusion (D aa * = < 
of a period is a superposition of the coherent 
states \a ) and | — a ). Such superposition 
states have interesting nonclassical properties 
and have been called Schrodinger cat states. 

Our stochastic treatment is designed to 
produce expectation values of operators for 
the oscillator. These can also be calculated 
analytically, but this is most easily performed 



I -Daa | — 2|/za 2 |) and hence is a good can- 
didate with which to test the ideas of Yuen 
and Tombesi. We should emphasize that the 
partial differential Eq. ([[£]) itself does not 
present any difficulties in spite of the nega- 
tive diffusion JTSJ. Indeed we can solve this 
equation directly to give the correct Q func- 
tion 



m 
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IV. ANALYTIC STOCHASTIC 

TREATMENT OF THE 
ANHARMONIC OSCILLATOR 

We can re-express the evolution of our Q 
function, 



given by Eq. (|i~9"D as an equiva- 
lent stochastic process using the method out- 
lined in Appendix [A]. A simple and natural 
choice is to set C aa * = = C a * a so that 

= v- 



C aa = \Ji2\ia and C a * a * = ^—i2fia*. The 
evolution equation for our Q function clearly 
displays negative diffusion and so we write 
our SDEs in terms of the variables a and a + . 
For the choices described above, our SDEs 
become 

a = —ifi2(a + a — l)a + £a (20) 
a + = ifi2{a + a - l)a + + (21) 

where £ and £ + are complex, white Gaussian 
noises with the stochastic averages 



(e(t)e(t')) s 
(t(t)m)s 



2ifi6(t - 1') 

-2ifi6(t - 1') 

0. (22) 



We will require averages over both the 
stochastic noise realizations and also over the 
initial quasi-probability distribution. The 
subscript 5* identifies the fact that we have 
carried out the stochastic average. The 
stochastic averages (22) do not fully deter- 
mine the forms of the noise terms. It is 
clear, however, that £ + (t) cannot be the com- 
plex conjugate of £(t). It has been suggested 
that the considerable freedom in choosing 
the forms of £(£) and £ + (£) can be used to 
suppress, although not completely remove, 
stochastic sampling errors in the analogous 
problem in the positive P representation. 
The analysis presented in this section is inde- 
pendent of this choice of Gaussian noise and 
hence the freedom to select the forms of £(t) 
and £ + (i) will not address the problem un- 
covered. 



We require the solution of Eqs. (|20|) and 
(pl|) with the initial conditions a(0) = (3 and 
a+(0) = p*. These mean that a+(0) = a*(0) 
and allow us to use the initial Q function to 
perform the average over the initial state. As 
already noted in Sect. [Qj, the form of the 
stochastic noise means that a + (t) will not 
take the value (a*(t)) in any given realiza- 
tion. The full quantum average will only be 
obtained by performing an average over the 
Qo function for the initial state. For the co- 
herent state |«o) this is 



-|/3-«o| 2 



7T 



(23) 



We denote the average obtained by integrat- 
ing over j3 by the subscript Qo- 

(F(a(t),a + (t))) Qo = [°° d 2 f3Q((3)F(a(t),a+(t) 



Quantum expectation values should be ob- 
tained on performing both the stochastic av- 
erage and the average over the this Qq func- 
tion. In particular, the mean value of a at 
time t will be 



(a(t)) = ((a(t)) s ) Qo . 



(24) 



We have not yet given a prescription for the 
order, if any, in which these averages must 
be performed. We will see that this question 
is of some importance for the solution of the 
SDEs. 

In this section we will calculate the ex- 
pectation value of the annihilation operator 
at time t by solving the SDEs (g^) and (0). 
We start by noticing that the combination 
a + a satisfies the equation 



— a + a — (£ + ^ + )a + a. 



(25) 



The formal solution to this equation is 

a + (t)a(t)=P*PeJ>' )+ t +it ' )]dt '. (26) 

This already suggests that the stochastic sim- 
ulation of this problem may run into difficul- 
ties. We expect the average obtained from 
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a + {t)a{t) will be (a(t)cfl(t)), which should 
take the constant value |«o| 2 +l- The solution 
(p6|), however, clearly shows that the stochas- 
tic noise will cause a + (t)a(t) to fluctuate 
away from its initial value in a single realiza- 
tion of the stochastic process. The average is 
constant but the corresponding variance in- 
creases in time. The presence of a complex 
driving force £ + means that a + (t)a(t) can 



acquire any complex value. Nevertheless we 
can proceed by inserting our solution (p6|) 
into our SDEs flSDp and ([21]). We find that 
the resulting equations are linear. In partic- 
ular, the equation for a(t) becomes 



a(t) 



the solution of which is 



a(t) = j3 exp 



i2»(\P\*e£M 1f W™*'-l)+W 



dt'y 



(27) 



Similar expressions have been given for the same model treated in the positive P represen- 
tation [0. The average of this quantity should be the expectation value of a(t). Let us 
start by performing the stochastic average. This can be achieved most readily by expanding 
the outer exponential in powers of |/3| 2 



Jo 



(a(t)) s = l3e i2 '*(ef W )dtf 

2/i 2 |/3| 4 f dt> C df'elo KM^Ml-eJT' KC-WOOI*' + . . . 



(28) 



Here we have made explicit use of the Gaussian nature of our stochastic noise in evaluating 
the averages of exponential functions of the noise. We can evaluate the average of each term 
in turn. The order zero and order one terms are 

(efi m *) s = e iia (29) 

-i2M 2 { [* dt'efo^+S + W ds eti^ s ' )ds )s = -*2/i|/3| 2 C dt'e 4 ^' e'^' 
Jo Jo 

= -\p\ 2 e lflt (e 2lflt -1). (30) 



It is straightforward to show that the 
stochastic average of the term of order \/3\ 2n 
is (-l) n \(3\ 2n e^\e 2i ^ -l) n /n\. It is tempting 
to re-sum the series in Eq. (EH) to give 



(a(t)) s = pe i3 ^ exp (|/5| 2 (1 - e i2 ^)) . (31) 

Let us see the consequences of this re- 
summing by completing our calculation of 
the expectation value of a(t) with the average 
over (3. This procedure leads to the expres- 
sion 



(Ht))s) Qo = - / d 2 (3e 

7T 



-|/3-qoI 2 ftJ3iJ,t 



Pe l6>lt exp 



Inspection of the integrand reveals a problem. 
It is clear that the integrand is unbounded 
(and the integral undefined) for times t such 
that cos(2/it) < 0. It is interesting to note 
that this includes the times, 7r/(2/i) and 
37r/(2/i), at which the anharmonic oscillator 
evolves into the Schrodinger cat states given 
in Eqs. ([ToD . The problem is that we have 
assumed that it is acceptable to perform the 
stochastic average before performing the av- 
erage over initial conditions. In fact this is 
not the case and we should perform the Qq 
, averaj^first. We can see this by evaluating 
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the average over the Qo function before re- for the exponential of a) a 0. We note that 
summing the series in our stochastic average this becomes the expression ([H]) obtained by 
given in Eq. (f28|). This gives the final average performing the stochastic average, if we iden- 
tify a and a) with f3 and (3* respectively. We 
have written Eq. (^) in antinormal order be- 
cause the Q function gives antinormally or- 



value 

(Ht)) Q o)s = 
oo 

= e «3M 



i2fit\n 



n=0 

oo 



E(i 



d 8 

— ^-/3 n+1 /3* ?1 e~' /3_a{ dered moments. If we use this expression to 
calculate the expectation value of a(t), for 



n=0 



E 



\a 



21 



1-- 



o(^ + l)! 



E 

1=0 



E(i 

n=l 



(n + 1) 



iO-ur initial coherent state, then we find 



+ l)\(n-l)\ 



a n oc 



(n + 1)! 
(n-l)W. 



(a(t)) =e 3 ^(a |S;exp((l 



e i2lxt )da) 



■\a ). (34) 



~ tfJjt a exp (Ja | 2 (e~ 



i2fit 



an, of course, evaluate this expectation 
which we recognize as the correct answer value by putting the operator into normal or- 
given in Eq. (|18|). Other moments can be dered form and using the fact that the coher- 
obtained in the same manner. ent states are right-eigenstates of the annihi- 

We can see the origin of the incorrect lation operator. Our aim, however, is to in- 
stochastic average given in Eq. ([31]) by con- vestigate the problems with the stochastic av- 
sidering the form of the annihilation operator erage associated with simulations designed to 
in the Heisenberg picture, Eq. (|i~6|) , which we reproduce antinormal ordered averages. We 
can also write in the form can work with the antinormal ordered form in 

e^f^aat; foo\ Eq. (|3~4D by expanding the exponential as a 
Taylor series and inserting the identity in the 
where : : denotes antinormal ordering and we form of an integral over the coherent states 
have used the antinormal ordering theorem | \p) [@]: 



a(t) 



3ifit ~ —i2fita^ a 

O CiC- 



e 3 ^aie (1 " 



(a(t)> 



(«ol E 



,n+l 



ra=0 



nl 



d 2 3 



7T 



\P){/3\a^\a ) 



oo l-\ i2fj,t\n 

n=0 



nl 



d 2 p 
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2n, 



l/5-ao| 2 



(35) 



Clearly it would be wrong to evaluate 
the summation before carrying out the in- 
tegral. Evaluating the integral first corre- 
sponds, in our stochastic treatment, to aver- 
aging over initial conditions before perform- 
ing the stochastic average and gives the cor- 
rect result. 

It is interesting to note that there is a 
strong similarity between the SDEs discussed 
here and those found for the anharmonic os- 
cillator in the positive P representation. In- 
deed, if we write equations for ae~ lflt and 
then we recover the equations dis- 



cre 



ifit 



cussed by Plimak et al [15|. An important 
difference, however, is that the diffusion for 



the positive P representation occurs with the 
opposite sign to that for the Q function. This 
means that the stochastic averages fl22|) have 
opposite signs when applied to the positive P 
representation. We can use the methods de- 
scribed in this section to obtain the expecta- 
tion value of a(t) in the positive P representa- 
tion. The stochastic average gives (a(t))s = 
f3e~ lflt exp (|/?| 2 (e~ ?2M * — 1)). Performing the 
average of this over a ^-function positive P 
distribution, peaked at /3 = a = (3 + *, gives 
the correct result ([Tj|). The positive P rep- 
resentation is associated with operator mo- 
ments in normal order and this seems to be 
the reason for the well-behaved form of the 



7 



stochastic averages for initial coherent states. 

V. STOCHASTIC SIMULATION OF 
THE ANHARMONIC OSCILLATOR 

In this section we present results of nu- 



merical simulations [19] of the stochastic pro- 
cess a(t) given in Eq. fl27|). Our simulations 
were performed using two discrete Gaussian 
processes r]i, rjf of the form 



Vi 



Af 



dtfffl 



Vi 



Af 



dt>H + (t f ) (36) 



where t = I At. In this way < tjitjii >= 
2ifxAt8iv and < r^r^ >= —2ifiAtSu>. We 
note that the relations do not completely 
specify the two independent complex white 
noises. As recently shown in the degree 
of freedom in the choice of the noise could 
be used to improve the results of the numeri- 
cal simulation by choosing the stochastic pro- 
cesses £ and £ + so as to inhibit (but not com- 
pletely suppress) the fast growth of a(t). In 
this paper, however, we have considered only 
the forms £ = \ / 2ifi(p and £ + = y/— 2ifi<p + 
with and <p + being real white noises. 

Each stochastic realization must start 
from a single point in phase space. For this 
reason, the analysis of the preceding section 
leads us to conclude that diverging trajec- 
tories, exploring large values of |a| are in- 
evitable. These divergences are responsible 
for the unbounded average obtained by per- 
forming the stochastic average before the av- 
erage over the initial Q distribution. Each 
of our simulations starts with at a point 
a(0) = a + *(0) = (3. Naturally, the aver- 
age over the initial Q distribution requires 
stochastic realizations for a range of values of 
P, weighted by the distribution (p3|). Consid- 
eration of a single value of (3, however, suf- 
fices to illustrate the divergences associated 
with individual trajectories. We observe, in 
each case, a divergence after some time. We 
can see the origin of these divergences in the 



SDEs (|20D and (gjj); the complex variables 
a and a + are not constrained to be complex 
conjugate quantities and so, in any given re- 
alization, the combination a + a can acquire 
an imaginary part. This leads to exponential 
growth of a or a + . The time at which this di- 
vergence appears varies between realizations 
and also depends on the initial conditions. 
In particular, the divergence appears earlier 
for larger values of \/3\ 2 . This is because of 
the exponential dependence of a(t) on |/3| 2 
as seen in Eq. (p7j). 

If we select a sufficiently small value of (3 
and perform an average over a large num- 
ber of trajectories then we find a result 
that is, for short times, in good agreement 
with the analytical average Eq. fl3Tf). In 
Fig. [I] we have plotted the time evolution 
of Re(a(t))s, obtained by considering 50000 
trajectories, starting from the initial condi- 
tion (3 = 0.001 + iO.l (diamonds line). For 
comparison the analytical expression for the 
stochastic average is represented by a contin- 
uous line. At very short times, we observe a 
near perfect agreement between the numeri- 
cal results and the analytical expression. At 
longer times, this agreement is lost because 
of the divergence induced by the independent 
stochastic noises. 

The trajectories start from a single point 
in phase space. This corresponds to select- 
ing, in each simulation, a 5 function phase- 
space probability distribution. Such a narrow 
distribution for the Qo distribution does not 
correspond to any physically allowed state 
@. Indeed, the evolution obtained from the 
Fokker-Planck equation for such an initial 
condition is highly singular. It is this behav- 
ior that is reflected in the divergent numerical 
simulations. Fig. ^| depicts the numerically 
obtained value of (a(t))s- We see that this 
average explores an extended region of the 
complex plane. The analytical average, Eq. 
(PI|), is represented by the small circle. 

The relationship between the time at 
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which trajectories diverge and the initial con- 
dition (/?) means that an ensemble of tra- 
jectories starting from a range of different 
initial conditions will rapidly produce diver- 
gences. For this reason the analytical result 
([UP cannot be reproduced numerically in any 
straightforward manner. 

VI. CONCLUSION 

In conclusion we have considered a pro- 
posal of Yuen and Tombesi |TU[ to give a 



stochastic representation of a Fokker-Planck 
equation with negative diffusion for the Q 
representation. We have shown that the cor- 
rect analytical moments for an anharmonic 
oscillator (associated with a non-linear 
process) can be obtained from the SDE's. 
These results, however, are highly sensitive to 
the order in which averages over the stochas- 
tic realizations and over the distribution of 
the initial conditions are performed. It is 
clear that more sophisticated techniques are 
required for stochastic simulation of the prob- 
lem. Recent work suggests that the effects 
of divergences can be significantly suppressed 
but not yet eliminated |15| , |20[| . 

The system studied in this paper is highly 
idealized. We could include the effects of loss 
and expect that these will improve the stabil- 
ity of the numerical results. Such an improve- 
ment has been noted in studies of the positive 
P 0. It is possible, however, that there are 
other interesting systems that are less sen- 
sitive to the noise and for these, stochastic 
simulations using the Yuen- Tombesi method 
may prove to be a useful technique. Possible 
systems for study in quantum optics include 
the Optical Parametric Oscillator and Second 
Harmonic Generation, that could be success- 
fully studied with this approach. Our prelim- 
inary studies suggest that there are regimes of 
operation, including the threshold, in which 
the probability for a divergent trajectory to 
occur is very small. In this case, numerical 



simulation does give stable results. We will 
return to this topic elsewhere. 
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APPENDIX A: 

In this Appendix we present a brief dis- 
cussion of the link between a given Fokker- 
Planck equation and an equivalent stochas- 
tic system (a more complete account can be 
found in @j2l|). As we have already noted, 
the Fokker-Planck equation does not corre- 
spond to a unique stochastic system and so 
it is natural to start with a stochastic system. 
Consider the pair of (Langevin) SDEs 

dot 

— = B a + C aa f a (t) + C aa *f a *(t) (Al) 
da* 

B a * + C a * a *f a *(t) + C a * a f a (t), (A2) 



dt 



with white Gaussian noise terms fi, defined 
to have zero average and second moments of 
the form 



</<(*)/i(* / )> = M(*-* / ) 



(A3) 



and the subscripts i,j denote a and a*. 
The formal solution of equations ( |A1|) and 
is: 



9 



a(t) = a(Q) + /* dt'Bjt') + 
Jo 



where the subscript po denotes an average 
over the initial probability distribution. The 
^^uantity p satisfies the conservation equation 

i} [dtp) = 0. (A7) 



9 9 ,■ s 



C aa dW(t') + f C aa *dW*(t') 
Jo 

a*(t) = a*(0) + r dt'B* a (t') + 
Jo 

J C a * a dW(t') + J C a * a *dW*(t'), (A5) The complete probability distribution is 

obtained by also averaging over the stochas- 
where we have introduced the Wiener pro- tic trajectories obtained with different noise 
cesses dW(t) = f a (t)dt and dW,(t) = realizations, denoted by the subscript [/(*)]: 
f a *(t)dt. 

In order to use these stochastic processes, 
we need to give a prescription for carrying 
out the stochastic integrals over the Wiener The time evolution for the distribution W 
processes. In this paper, we choose the can be obtained using the continuity equa- 
Stratonovich interpretation of the stochastic tion and gives [[4|,|22||: 
integral in which 

Emu) - w(u .) ■ | w - ~i B ^ - w + \ E ^ 



W(a,a*,t) = (p(a,a*,t; [f(t)\ )) [/(t)] . (A8) 



g(a,a*)dW(t') 



1. 1 1 



g | a (tj) + a (tj-i) a (U) + a (^~i)w here the subscripts k again denote a 
\ 2 2 a *. if compare this form of the Fokker- 

The reason of this choice, instead of the Ito Planck ec i uation with the Eq. © then we 
interpretation, is that we will be constructing obtain the correspondences: 



analytical averages over the stochastic pro- 
cess and the Stratonovich formalism allows 
us to use the familiar rules of calculus. 

From the Langevin equations it is pos- 
sible to obtain a unique Fokker-Planck 
equation for the probability distribution 
W(a(t), a*(t), t). If we consider the tra- 
jectory obtained in a single realization 
of the stochastic process / = (fa, fa*) 
and start from the initial value 5(0) = 
(a(0), a*(0)), then the solution at time t 
is completely determined and the proba- 
bility distribution for it is the 5-function 



Do 



D, 

A a = B a 
( d 



r 2 + r 2 
= r 2 + r 2 

C a aC a *a C a *a*C a 



da 



Cn 



1 

2 

a 



d_ 

da 



C aa I C aa ~t~ 



' d_ 

da* 



_d_ 

da" 



Co 



(A9) 
(A10) 
(All) 



(A12) 



A, 



Bo 



1 

+ 2 



^~>a*a 

. da 



C 1 aa 



'A 

da* 



Ca*a ) C a *a~\~ 



d_ 

da 



Ca*a* l C a 



+ 



'_d_ 

da* 



(A13) 



5 (a - a(t; 5(0); [/(*)])) 5 (a* - a*(t; 5(0); [f(t)\)). 
Considering a set of initial conditions 5 , dis- Note that these equations do not give unique 
tributed according some initial distribution forms for the B and C functions. This is a 
Po = p(a(0);0), we can obtain the shape of consequence of the lack of a unique stochas- 
the distribution at time t: tic representation for a given Fokker-Planck 



p(a,t; [f(t)}) = 
(5 (a - a(t- 5(0); \f(t)))) 5 (a* - a*(t; 5(0); [f(i 



equapon. 



our stochastic variables a and a* are 
Picomplex conjugate quantities, then it 
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follows from equations ( |A1|) and (|A2|) that 



C n * n and C* 



Crusty* . 



These conditions necessarily imply positive 
diffusion as, from ( |A9| ) to (|A11|) , D aa * > 
\D aa \. It follows that the stochastic variables 
cannot be complex conjugate quantities when 
we have negative diffusion. In order to avoid 
possible confusion, we replace the stochastic 
variable a*(t) by a + (t) whenever there is neg- 
ative diffusion. 
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FIG. 1. Time evolution of the stochastic a- 
verage Re(a(t))s, using 50000 trajectories and 
starting from (3 = 0.001 + £0.1. The diamonds 
represent numerical values. The continuous line 
represents the analytical result for (a(t))s- The 
dashed line represents {{a(t))Q )s for the initial 
coherent state |«o > with oq = 0.001 + £0.1. 
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FIG. 2. Phase space plot of the numerical 
average {a(t))s, the real part of which is shown 
in Fig. |]. The points represent numerical val- 
ues at different times. The circle represents the 
analytical result for (a(t))s- 
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